{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6f6f2508",
   "metadata": {},
   "source": [
    "# Robustness Check for PSM\n",
    "\n",
    "The results of the propensity score matching estimation method can be unstable (Kind and Nielsen,2019). To verify the degree of instability of the results, we have added the following functions to help everyone identify this, including: \n",
    "* Calculate the total number of observations remained, you should avoid discard too many observations.\n",
    "* Calculate the most repeated times of a control observations, if a control observation appears too many times, something might be wrong for the propensity score estimation model.\n",
    "* Sensitivity test, following Rosenbaum (2005) \"Sensitivity analysis in observational studies\". You can start with fitting a list of integers starting from 1, such as [1,2,3,...] and see from which gamma, the \"p-val upper bound\" start to fall below 0.05.\n",
    "* Placebo test, verifies if the ate disappears when the treatment is randomly assigned.\n",
    "\n",
    "\n",
    "### Data\n",
    "\n",
    "We have following types of observations:\n",
    "* Covariates, which we will denote with `X`\n",
    "* Treatment, which we will denote with `T`\n",
    "* Responses, which we will denote with `Y`\n",
    "\n",
    "Requirement is that `T` is a binary varible which contain only 0/1 values. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7c24ac12",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['/Library/Frameworks/Python.framework/Versions/3.12/lib/python312.zip', '/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12', '/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/lib-dynload', '', '/Users/bytedance/Library/Python/3.12/lib/python/site-packages', '/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages']\n",
      "current version is:  0.0.5\n",
      "/usr/local/bin/python3.12\n",
      "3.12.2 (v3.12.2:6abddd9f6a, Feb  6 2024, 17:02:06) [Clang 13.0.0 (clang-1300.0.29.30)]\n",
      "sys.version_info(major=3, minor=12, micro=2, releaselevel='final', serial=0)\n"
     ]
    }
   ],
   "source": [
    "# import sys\n",
    "# print(sys.path)\n",
    "# sys.path.append('/Users/bytedance/PycharmProjects/github/CausalMatch')\n",
    "import causalmatch as causalmatch\n",
    "from causalmatch import matching,gen_test_data\n",
    "print('current version is: ',causalmatch.__version__)\n",
    "import sys\n",
    "print(sys.executable)\n",
    "print(sys.version)\n",
    "print(sys.version_info)\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "\n",
    "df, rand_continuous, rand_true_param, param_te , rand_treatment, rand_error = gen_test_data(n = 10000, c_ratio=0.5)\n",
    "X = ['c_1', 'c_2', 'c_3', 'd_1', 'gender']\n",
    "y = ['y', 'y2']\n",
    "id = 'user_id'\n",
    "\n",
    "# treatment variable has to be a 0/1 dummy variable\n",
    "# if is string, please convert to a 0/1 int input\n",
    "T = 'treatment'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "55b58610",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -------------------------------------------------------- #\n",
    "# STEP 1: initialize matching object\n",
    "match_obj = matching(data = df,\n",
    "                     T = T,\n",
    "                     y = ['y'],\n",
    "                     X = X,\n",
    "                     id = id)\n",
    "\n",
    "# -------------------------------------------------------- #\n",
    "# STEP 2: propensity score matching without trim anything\n",
    "match_obj.psm(n_neighbors = 1, \n",
    "              model = LogisticRegression(), \n",
    "              trim_percentage = 0, \n",
    "              caliper = 1) \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3623eb97",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Wilcoxon-statistic</th>\n",
       "      <th>gamma</th>\n",
       "      <th>stat upper bound</th>\n",
       "      <th>stat_lower_bound</th>\n",
       "      <th>z-score upper bound</th>\n",
       "      <th>z-score lower bound</th>\n",
       "      <th>p-val upper bound</th>\n",
       "      <th>p-val lower bound</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>7670016.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>6496126.50</td>\n",
       "      <td>6496126.50</td>\n",
       "      <td>11.173314</td>\n",
       "      <td>11.173314</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>y</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>7670016.0</td>\n",
       "      <td>1.5</td>\n",
       "      <td>7795351.80</td>\n",
       "      <td>5196901.20</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>y</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7670016.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>8661502.00</td>\n",
       "      <td>4330751.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>y</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7670016.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>9280180.71</td>\n",
       "      <td>3712072.29</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>y</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>7670016.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>9744189.75</td>\n",
       "      <td>3248063.25</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>y</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Wilcoxon-statistic  gamma  stat upper bound  stat_lower_bound  \\\n",
       "0           7670016.0    1.0        6496126.50        6496126.50   \n",
       "1           7670016.0    1.5        7795351.80        5196901.20   \n",
       "2           7670016.0    2.0        8661502.00        4330751.00   \n",
       "3           7670016.0    2.5        9280180.71        3712072.29   \n",
       "4           7670016.0    3.0        9744189.75        3248063.25   \n",
       "\n",
       "   z-score upper bound  z-score lower bound  p-val upper bound  \\\n",
       "0            11.173314            11.173314                0.0   \n",
       "1                  NaN                  NaN                NaN   \n",
       "2                  NaN                  NaN                NaN   \n",
       "3                  NaN                  NaN                NaN   \n",
       "4                  NaN                  NaN                NaN   \n",
       "\n",
       "   p-val lower bound  y  \n",
       "0                0.0  y  \n",
       "1                NaN  y  \n",
       "2                NaN  y  \n",
       "3                NaN  y  \n",
       "4                NaN  y  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# -------------------------------------------------------- #\n",
    "# STEP 3: Sensitivity test, following Rosenbaum (2005) \"Sensitivity analysis in observational studies\".\n",
    "#         You can start with fitting a list of integers starting from 1, such as [1,2,3,...] and\n",
    "#         see from which gamma, the \"p-val upper bound\" start to fall below 0.05\n",
    "#         !!!!! only works for PSM  !!!!!\n",
    "\n",
    "gamma = [1, 1.5, 2, 2.5, 3]\n",
    "match_obj.sensitivity_test(gamma = gamma)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a224c9a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "empirical critical value:  -0.31556721533936816 0.2844152580884073\n",
      "if pass placebo test:  Pass\n"
     ]
    }
   ],
   "source": [
    "# -------------------------------------------------------- #\n",
    "# STEP 4: placebo treatment, verifies the effect disappears when the Treatment is replaced with a placebo.\n",
    "#         Returen a numpy array to store ATE estimate for each bootstrap sample with length b X 1, default is 100.\n",
    "ate_array = match_obj.placebo_treatment()\n",
    "print('empirical critical value: ',\n",
    "      match_obj.critical_value_lb, \n",
    "      match_obj.critical_value_ub)\n",
    "\n",
    "\n",
    "b_pass_placebo_test = 'Fail' if (match_obj.critical_value_lb <= match_obj.ate().iloc[0]['ate'] and match_obj.critical_value_ub >= match_obj.ate().iloc[0]['ate']) else 'Pass'\n",
    "print('if pass placebo test: ',b_pass_placebo_test)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4f0e1709",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             Robustness Check Output Table for Dep. Variable y              \n",
      "============================================================================\n",
      "                                                            coef    P>|t|   \n",
      "  Average Treatment Effect:                                 0.52     0.00%  \n",
      "                                                                            \n",
      "  1. Total % of obs remained:                            101.94%      Pass  \n",
      "   -- Treated % obs remained:                            100.00%      Pass  \n",
      "   -- Control % obs remained:                            103.96%      Pass  \n",
      "      -- The most repeated times of a control obs:             9         -  \n",
      "  2. Sensitivity test result, Gamma statistics               1.0      Fail  \n",
      "  3. Placebo Test Result conf. interval              -0.32,0.28,      Pass  \n"
     ]
    }
   ],
   "source": [
    "# -------------------------------------------------------- #\n",
    "# STEP 5: Robustness check \n",
    "#  -- \"1. Total % of obs remained: \" is used to check if you discard too many observations or not\n",
    "#         in this example, \n",
    "#  --    \"Control % obs remained: \" is greater than 100% because we set \"trim_percentage = 0\"\n",
    "#         and \"caliper = 1\", which means we keep all observations\n",
    "#  --    \"The most repeated times of a control obs:\" means for an observation in control group, how many\n",
    "#        times it has been used to match to different treated observation. If an observation has been\n",
    "#        used too many times, something might be off with your propensity score model or your dataset.\n",
    "\n",
    "match_obj.robust_check(gamma)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
